adams_bashforth_coefficients(1, 1)
adams_bashforth_coefficients(2, 2)
adams_bashforth_coefficients(3, 3)
adams_bashforth_coefficients(4, 4)

function adams_bashforth_coefficients(s, p)
    % 定义符号变量 b0, b1, ..., b_{s-1}
    syms b [1 s]  % 创建 b0, b1, b2, ..., b{s-1}（MATLAB 下标从 1 开始）
    syms bs
    eqn = bs == 1;
    % 设置阶条件方程
    eqns = sym(zeros(1, p));  % 需要 p 个方程
    for m = 1:p
        j = 0:s-1;  % j = 0, 1, ..., s-1
        if m == 1
            eqns(1) = sum(b) == 1;  % C_1 条件: b0 + b1 + ... + b_{s-1} = 1
        else
            % 高阶条件: sum(j.^(m-1) * b_j) = [s^m - (s-1)^m] / m
            eqns(m) = sum(j.^(m-1) .* b) == (s^m - (s-1)^m)/m;
        end
    end
    
    % 求解方程组
    sol = solve(eqns, b);
    sols = solve(eqn, bs);
    % 提取并显示系数（转换为分数形式）
    fprintf('\nAdams-Bashforth 系数 (s=%d, p=%d):\n', s, p);
    
    % Handle different solution formats
    if s == 1
        % For s=1, sol is returned directly as a value
        fprintf('b0 = %s\n', char(rats(double(sol))));
        fprintf('b1 = %s\n', char(rats(double(sols))));
    else
        % For s>1, sol is a structure
        for k = 1:s
            bk = sol.(sprintf('b%d', k));
            fprintf('b%d = %s\n', k-1, char(rats(double(bk))));
        end
        fprintf('b%d = %s\n', s, char(rats(double(sols))));
    end
end